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ABSTRACT 


A semi-analytical method for determining the strain energy release rate due to a prescribed 
interface crack in an adhesively-bonded, single-lap composite joint subjected to axial tension is 
presented. The field equations in tenns of displacements within the joint are formulated by using 
first-order shear defonnable, laminated plate theory together with kinematic relations and force 
equilibrium conditions. The stress distributions for the adherends and adhesive are determined 
after the appropriate boundary and loading conditions are applied and the equations for the field 
displacements are solved. Based on the adhesive stress distributions, the forces at the crack tip 
are obtained and the strain energy release rate of the crack is determined by using the virtual 
crack closure technique (VCCT). Additionally, the test specimen geometry from both the ASTM 
D3165 and D1002 test standards are utilized during the derivation of the field equations in order 
to correlate analytical models with future test results. The system of second-order differential 
field equations is solved to provide the adherend and adhesive stress response using the symbolic 
computation tool, Maple 9. Finite element analyses using /-integral as well as VCCT were 
performed to verify the developed analytical model. The finite element analyses were conducted 
using the commercial finite element analysis software ABAQUS™. The results determined 
using the analytical method correlated well with the results from the finite element analyses. 

KEY WORDS: adhesively-bonded joint, laminated plate theory, composite joint, strain energy 
release rate, virtual crack closure technique (VCCT), fracture mechanics 
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applied tensile force per unit width, N/m 

x-directional displacement, m 

mid-plane x-directional displacement, m 

bending slope, radians 

z-directional displacement, m 

normal stress resultants per unit width, kN/m 

bending moment per unit width, kN 

transverse shear stress resultant per unit width, kN/m 

in-plane modulus per unit width, kN/m 

transverse modulus per unit width, kN/m 

coupling modulus per unit width, kN 

flexural modulus per unit width, kNm 

shear correction factor 

adhesive shear stress, kPa 

adhesive peel stress, kPa 

adherend thickness, m 

adhesive thickness, m 

adhesive normal stress, kPa 

adhesive strain 

Young’s modulus of adhesive, kPa 
shear modulus of adhesive, kPa 
Poisson’s ratio of adhesive 
length of prescribed crack, m 
length of virtual crack extension, m 
overlap length before crack is initiated, m 
notch length of ASTM D3 165 specimen, m 
crack tip forces 

work required to close crack propagation b per unit joint width, N 

total strain energy release rate, N/m 

crack tip forces from finite element model, N 

./-integral value 
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INTRODUCTION 


Advanced composite materials have been widely used due to their high strength-to- 
weight ratio and excellent corrosion resistance. In many applications, bolted joints have been 
replaced by adhesively-bonded joints because of the weight penalty and corrosion problems 
associated with bolted joints. However, the geometric discontinuity present at the ends of an 
adhesively-bonded joint results in peak shear and normal (peel) stresses in the adhesive layer that 
typically occur near the adhesive to adherend interface. The stress concentration caused by the 
discontinuity often results in local yielding and further develops into crack initiation. In a 
research study conducted by Yang et al. [1] on adhesively-bonded joints, they concluded that the 
fracture mechanics approach would be an effective method for predicting the load carrying 
capacity of a bonded joint. 

Earlier studies of adhesively-bonded joints can be found in the extensive reviews given 
by Kutscha [2], Kutscha and Hofer [3], Matthews et al. [4], and Vinson [5]. Yang and Pang [6] 
derived an analytical model that provided the stress distributions of adhesively-bonded single-lap 
composite joints subjected to axial tension. Huang et al. [7] and Yang et al. [8] also derived an 
elastic-plastic model for adhesively-bonded single-lap composite joints. Important capabilities 
included in their approaches were the asymmetry of the adherend laminates as well as the effects 
due to transverse shear deformation. 

An existing crack is usually assumed in a joint when conducting a fracture analysis. 
Krueger [9] described the virtual crack closure technique (VCCT), including its history, 
approach, and applications, in conjunction with finite element analysis in a report published in 
2002. Davidson et al. [10-13] published a series of papers that employed the classical plate 
theory version of the VCCT to predict the strain energy release rate of mixed-mode delamination 
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in composite laminates. A crack-tip force method was derived by Park and Sankar [14] to 
compute the strain energy release rate in delaminated beams and plates. Kim et al. [15] proposed 
a simplified method for detennining the strain energy release rate of free edge delamination in 
composites using the classical laminated plated theory. Finite element methods play a significant 
role in structural analysis and have been widely used to study the adhesively-bonded composite 
joint. Wang et al. [16] applied the VCCT to calculate the strain energy release rate of cracked 
composite panels with nonlinear deformation. Wei et al. [17] presented an improved VCCT to 
determine the energy release rate using a three-step analysis. Yang et al. [1] developed finite 
element models using the finite element software, ABAQUS™, to estimate the /-integral of an 
adhesively -bonded joint with a crack. 

Although finite element analysis methods are capable of solving problems with varying 
material types and complicated geometrical configurations, analytical solutions offer 
performance and solution advantageous, especially when performing parametric analyses and 
optimization. The objective of the present paper is to describe a semi-analytical fracture 
mechanics method that can be used to determine the strain energy release rate due to a prescribed 
crack in an adhesively-bonded single-lap composite joint. The prescribed crack is assumed to be 
present at the location of peak stress in the overlap region, usually at the corners of the adhesive 
adjacent to the continuous adherend. Additionally, the strain energy release rate calculated using 
this method can be applied to configurations with adhesive or adhesive/adherend interfaces as a 
failure criterion given the same mixed-mode critical strain energy release rate exists in the 
subject joint to account for the adhesive cohesion or adhesion failure. Linear elastic material 
properties as well as small displacements are assumed for both the adhesive and adherends in 
order to make the analytical approach feasible. Due to the fact that the geometric nonlinearity 
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during the joint deformation is not considered, the analytical model is more accurate when 
analyzing joints with stiffer and/or thicker adherends. 

In the remaining sections, the displacement fields and stress distributions of adhesively- 
bonded single-lap composite joints with a prescribed interface crack under tension are 
determined analytically using laminated anisotropic plate theory. After the virtual crack 
extension is applied, the crack tip forces are calculated based on the adhesive stress distributions. 
The strain energy release rate is then calculated using the VCCT. ASTM D3165 [18] and ASTM 
D1002 [19] specimen geometries are used in the model derivations. The semi-analytical 
solutions are determined using the symbolic computation tool Maple 9 [20]. Results from the 
analytical model are verified by finite element analysis using ABAQUS™ 6.3 [21, 22], 

ANALYTICAL DEVELOPMENT OF FIELD EQUATIONS FOR THE 
CALCULATION OF ENERGY RELEASE RATES 

The details of the semi-analytical method for determining the strain energy release rate of 
adhesively-bonded single-lap composite joints with a prescribed interface crack are now given. 
Based on the laminated plate theory version of Irwin’s virtual crack method [9], the strain energy 
release rate is derived in terms of the forces and moment at the crack tip, N c , Q c , and M c . These 
forces and moments at the crack tip are determined from the linear elastic shear and peel stress 
distributions in the adhesive. Therefore, a description of the stress state in the pre- and post- 
cracked specimen geometry is required before an estimate of the strain energy release rate can be 
obtained. 

ASTM D3165 Stress Model 

A summary of the basic methodology used by the authors to derive the equations for 
determining the required stress and displacement fields in an adhesively-bonded joint are 
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presented in the present section. The details of the method used to determine the displacement, 
strain, and stress fields in an adhesively bonded joint are available in the previous report by Yang 
and Pang [6]. Shown in Figure 1 is an adhesively-bonded single-lap joint with the standard 
geometry of an ASTM D3165 specimen and an applied tensile load P per unit width. The 
specimen geometry is divided into five regions to aid in the derivation of model equations. 
Region 3 is the bonded joint overlap area where the applied mechanical loads are transferred 
from one adherend to the other, and is also the area for which the joint strength is typically based. 

The generalized field equations for the adherends and adhesive are the same for all the 
three regions with an adhesive layer. The behavior of the adherends is described by the 
laminated anisotropic plate theory and the adhesive is assumed to behave as an elastic, isotropic 
material. The displacement fields for the upper and lower adherends can be written as follows: 


u = u° u (x) + z u If/ U (x) (1) 

u L = u oL {x) + z L iff L {x) (2) 

w = w (x) (3) 

w L = w L {x) (4) 


where displacements in the axial (x) direction are given by u, displacements in the transverse (z) 
direction are given by w, the superscripts U and L denote the upper and lower adherends, 
respectively, superscript o represents the mid-plane displacement, and i// is the corresponding 
bending slope. The normal stress resultant per unit width, N x , bending moment per unit width, 
M y , and transverse shear stress resultant per unit width, Q z , are expressed as 
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Figure 1 . ASTM D3165 specimen configuration and coordinate systems. 
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where z a and q a are the shear and peel stresses of the adhesive, respectively, and h lJ is the 
thickness of the upper adherend. Three equilibrium equations can be obtained for the lower 
adherend in a similar fashion to show: 


dN x _„ 

dx 


( 11 ) 
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Figure 2. Free body diagram and sign conventions. 

Using the kinematics of the adherends and assuming a perfect bond between the adhesive 
and the adherend surfaces, the adhesive strains are related to the bottom surface of the upper 
adherend and the top surface of the lower adherend. In terms of the displacement field of the 
two adherends, the adhesive strains in the traditional sign convention with x, v, z subscripts can 
be written as: 


r. 


n 


( oU oL\ [ h L U , h L L 

\u -U j-| — v ^ 


dw _j_ dw 


U\ 


dx dx 


( 14 ) 
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where ;/ is the adhesive thickness, and h l and /? £ are the thickness of the upper and lower 
adherends, respectively. Assuming a condition of plane-strain, the adhesive stresses can be 
obtained as: 



(17) 

VE “ ( + ) 
cry /. y, ~ £:) 

(l+vjl-2vj 

(18) 

£ 

q ° =C 7 ~n+ vi + 

(\+v a )(\-2v a ) 

(19) 

Ta = - G a r xz 

(20) 


As previously noted, r a and q a are the shear and peel stresses of the adhesive, respectively, and 
their sign conventions are shown in Figure 2. E a is the Young’s modulus, v a is the Poisson’s 
ratio, and G a is the shear modulus of the adhesive. 

Substituting the stress resultants in Equations (5) - (7) for the upper and lower adherends 
along with the adhesive stresses in Equations (19) and (20) into Equations (8) - (13), six coupled 
second-order ordinary differential equations with six unknowns u oU , t// u , w u , u° L , i//, and w L are 
obtained for each of Regions 1,3, and 5. In Region 2, only the upper adherend is present 
without any adhesive stresses. Therefore, the force equilibrium conditions in Region 2 are 
similar to those in Equations (8) - (10) but without the adhesive stresses, z a and q a , 

^ = 0 ( 21 ) 

dx 
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( 22 ) 


dMy _ U 

dx 



= 0 


( 23 ) 


Substituting the stress resultants in Equations (5) - (7) of the upper adherend into Equations (21) 
- (23), three coupled second-order ordinary differential equations with three variables u oU , i//\ 
and w u are obtained for Region 2. Applying the similar approach to Region 4 yields another 
three coupled second-order ordinary differential equations with three variables u oL , i//, and w L . 

The overall system of governing equations, including all five regions, contains twenty- 
four second-order ordinary differential equations with twenty-four unknown variables. A total of 
forty-eight boundary conditions are obtained at the two ends of each Region based on either 
continuity or applied force conditions. The symbolic solver Maple 9 was used to solve the 
system of equations and obtain the displacement, strain, and stress fields. 


Strain Energy Release Rate Calculation 

Once the stress, strain, and displacement fields are known in the adhesively-bonded 
single-lap composite joints, the virtual crack closure technique (VCCT) is applied to estimate the 
strain energy release rate of the joint with a prescribed interface crack. According to linear 
elastic fracture mechanics, the energy released from the propagation of a small crack is 
equivalent to the work needed to close that small crack propagation [16]. 

When a crack is initiated, it usually starts at a location of high stress concentration. For 
the ASTM D3165 specimen configuration shown in Figure 1, the critical areas are located at the 
upper left and lower right corners of the adhesive layer in Region 3. In deriving an expression 
for the strain energy release rate, a joint is assumed to have an overlap length l Q , a notch size /„, 
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and a crack of length a that is located at the lower right adhesive/adherend interface. Due to the 
anti-symmetry of the joint, only the right half of the joint, as shown in Figure 3, is used for the 
strain energy release rate calculation. The displacement of the crack tip C after a load is applied 
can be determined using the stress model previously described and a Region 3 length of L 2 = l 0 - 
2 a and Regions 2 and 4 lengths of L 2 = L 4 =/„ + a. In the stress fonnulation, two boundary 
conditions used for solving the field equations are that the displacements in both the x and z 
directions of the left end of the joint (jci = 0) are zero. Because only the right-half of the joint is 
used in the strain energy release rate calculation, the center point O of the joint, as shown in 
Figure 3, is used as the reference point to avoid double-counting the effects due to the crack. 
Both the displacements in x and z directions, a and w, in the strain energy release rate calculation 
are relative displacements to point O. However, the absolute bending slope, (//, is used directly 
for the strain energy release rate calculation. According to the displacement fields for the upper 
and the lower adherends, the relative displacements of the crack tip C to the reference point O 
are obtained using the stress model and a joint with a central overlap length L 2 = l 0 - 2 a and 
notch lengths L 2 =L 4 = /„ + a, i.e. 


^ n "1 



IJ 2 


Figure 3. A joint of central overlap length l with a prescribed interface crack of length a. 
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U C,0~ U \x 3 =l„-2a + 2 ^ I x 3 =l 0 —2a U 0 


YC-Y \x,=l„-2a 


W C,Q= W \x 3 =L-2a~ W 0 


where uqo and wc.o are the relative displacements of crack tip C to the reference point O in the x 
and z-directions, respectively, xj, = l 0 - 2 a denotes the location of the crack tip C. The 
displacements, uo and wo, of the center reference point O, are determined from the average 
displacements of the bottom surface of the upper adherend and the top surface of the lower 
adherend at vy = (l 0 -2a)/2 as, 


u 0 =~\ u 


1 ( h u h L 

2 U \x 3 =(l„-2a)/2 \x } =(l 0 -2a)/2 +U \x 3 =(I 0 -2a)/2 \x 3 =(l 0 -2a)/2 


_ 1 U | . I | 

W 0~~^\ W \x } =(lo-2a)H +W \x 3 =(l„-2a)l2 j 


If the crack propagates a small length b (virtual crack), the crack tip moves to point C\ and the 
previous crack tip C separates into points A and B , as shown in Figure 4. 


Upper Adherend 


Adhesive 


1 1 C’ s 3 

'■ in t I rWKtpy AAlnpypyiA 



h Lower Adherend Crack Tip 



b 1 a 


Figure 4. Virtual crack extension of length b. 
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Before the crack growth, the adhesive between C’ and C adheres to the lower adherend 
where adhesive shear and peel stresses exist at the interface as shown in the left of Figure 5. If 
the crack propagates a small length b, the adhesive shear and peel stresses vanish between C’ and 
C. The crack tip forces corresponding to a small crack propagation b can be assumed as the 
equivalent forces, Nc, M c , and Q c , to the shear and peel stresses between C’ and C, as shown in 
the right of Figure 5, and can be calculated as 



Figure 5. Equivalent forces at the crack tip. 

N C = j/ f . T a dX i ( 29 ) 

Jl 0 -2a-b 

M c = q “ ( /o ~~ 2a " x ’ ) Jx 3 ( 30 ) 

Q C = “I/' ^ ,fa dX 3 (31) 

where the adhesive shear stress z a and peel stress q a are obtained from the stress model with the 
overlap length Lj = l a - 2a. Note that only the upper adherend and adhesive layer with point A 
are shown in Figure 5 and that the positive directions of Nc, Me, and Qc are defined to be 
consistent with the positive directions of displacements u, i //, and w, respectively. The same 
approach is applied to point B to obtain the crack tip forces at the lower adherend. It can be seen 
that the crack tip forces at point A have the same magnitude but opposite directions as those at 
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point B because the adhesive stresses acting on the upper adherend are in the opposite directions 
as those acting on the lower adherend. 

After the crack propagates a small length b, the previous crack tip C separates into two 
points A and B , as shown in Figure 4. To determine the relative displacements of point A and B 
to the reference point O, a subsequent stress analysis is performed using a joint with a central 
overlap length = l 0 - 2a - 2b, which simulates the overlap up to the new crack tip C\ 
including the two notch lengths L 2 = L 4 = l n + a + b. Extension of the upper adherend from C ’ to 
A with a free end at A generates the displacements at point A. Therefore, the relative 
displacements of point A to the reference point O are 




( 33 ) 

( 34 ) 

( 35 ) 

( 36 ) 


Wr.o = w " U=b ~ w o ( 37 ) 

The subscripts A and B denote points A and B. All the other superscripts and subscripts are the 
same as previously noted. Because of the different joint configuration in the stress model before 
and after the crack propagates, the displacements of the reference point O used in Equations (29) 
- (34) are obtained again from a stress analysis with a central overlap length L 3 = l Q - 2a - 2b as 
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\x i =(l a -2a-2b)/2 


\x } =U„-2a-2b)/2 , 


The central point displacements, u 0 and w 0 , as shown in Equations (27) and (28) are before the 
crack propagates and the central point displacements, uj and w„\ as in Equations (38) and (39) 
are after the crack propagates. In order to close the virtual crack propagation, of length b, the 
crack tip forces are applied on points A and B to move them back to the location of the original 
crack tip C. Therefore, the total work required to close the small crack propagation b is 


W = | [N c (u co -u ao )+M c ( y/ c -y/ A )+ Q c (w co - w A0 )] 

~ T [^c( W C,0 _ li B,0 )"*" ^ C^\V C — W B ) "*“ 0c( W C,O ~ 


After simplification, the work becomes 


W— ^N c {u bo u ao )+M c (i// b Ya ) + Qc( w b,o w a,o )] 


For joints with a unit width, the strain energy release rate is defined as the derivative of energy 
released from the crack propagation with respect to the length of the crack propagation 



( 42 ) 


where U is the strain energy stored in the body. Based on the VCCT, the total energy released 
from the crack propagation is equivalent to the work needed to close the same crack propagation. 
With a virtual crack propagation of length b, the total strain energy release rate GY, which is a 
summation of the Mode I strain energy release rate G/ and the Mode II strain energy release rate 
Gu, can be calculated as 


16 



G T =G I +G II =^-= J ^- = \N c (u b - u A ) + M c {y/ B -y/ A )+ Q c (w A - w B )] (43) 

da b 2b 

FINITE ELEMENT MODEL DESCRIPTIONS USING 
VCCT AND /-INTEGRAL METHODS 

Finite Element Models using VCCT 

The strain energy release rate due to a small increase in crack length is equivalent to the 
energy required to close that small crack increment. Therefore, the strain energy release rate can 
be computed by finite element models using the VCCT [15], As shown in Figure 6, the tip of a 
crack with an original length a is located at C. Assuming a virtual crack propagation of length b, 
the new crack tip becomes C’ and the original crack tip becomes two separate nodes/ and g. If 
nodes/ and g are restrained at the original crack tip location, this virtual crack with length b is 
closed and the work to close this virtual crack can be calculated by multiplying the reaction 
forces at nodes / and g by the relative displacements of these two separate nodes to the original 
crack tip C. 

f 

'3 


a b 

Figure 6. Demonstration of finite element method with VCCT. 

The virtual crack closure technique originally developed by Rybicki and Kanninen [24] as well 
as Raju [25] and further developed by Wei [17] is applied in the current study. For 2-D 
conditions, the following steps are required: 
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1. Build the finite element model with an original crack of length a and detennine the 
displacements of crack tip C, u c and w c in the x and z directions, respectively. 

2. Propagate the crack with a small length b (usually one element size); the original crack 
tip C becomes two separate nodes. Record the displacements of both nodes, Uf, Wf, u g and 
w g . 

3. Constrain the two separate nodes so that they have the same displacements as the original 
crack tip C and obtain the reaction forces F x f, F y f, F xg , and F yg . 

4. The work needed to close the virtual crack is 

w = \ i F xf (« c - u f) + F A w c~ w f)Y \ K Oc - u g ) + F yg O c - w g )] (43) 

The total strain energy release rate is then obtained by 

W 

0, = ^ (44) 

b 

In the present study, two-dimensional 4-node linear plane-strain quadrilateral elements 
are utilized in the finite element model for VCCT application. There are 12 elements through the 
adherend thickness corresponding to 12 plies of T300/5208 with orientation and sequence of 
[0 3 /90 3 ]s and 6 elements are used in the finite element mesh through adhesive thickness. 

Finite Element Model Description using the /-Integral 

Finite element models with the /-integral calculation were constructed using ABAQUS™ 
[21, 22] to verify the present analytical model. The /-integral is usually used in quasi-static 
fracture analysis to characterize the energy release associated with crack propagation. It is 
equivalent to the strain energy release rate if the material response is linear elastic. Considering 
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an arbitrary counter-clockwise path (F) around the crack tip, as illustrated in Figure 7, the J- 


integral is defined as 


J = J wdy-Y J T i~^ Lds 


where w is the strain energy density, uj, u 2 , and u 3 are the components of the displacement vector, 
ds is the incremental length along the contour F and 77, T 2 , and 77 are components of the 
traction vector. The traction is a stress vector normal to the contour. In other words, 77, 77, and 
77 are the nonnal stresses acting at the boundary if a free-body diagram on the material inside of 


the contour is constructed. 



Figure 7. An arbitrary contour around crack tip. 

Several contour integral evaluations are possible at each location along the crack front. 
In a finite element model each evaluation can be thought of as the virtual motion of a block of 
material surrounding the crack tip. Each such block is defined by contours and each contour is a 
ring of elements completely surrounding the crack tip or crack front from one crack face to the 
opposite crack face. These rings of elements are defined recursively to surround all previous 
contours. ABAQUS/Standard automatically finds the elements that form each ring from the 
node sets given as the crack-tip or crack-front definition. Each contour provides an evaluation of 
the contour integral [22], 
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Theoretically, the ./-integral should be independent of the domain used, but ./-integral 
estimated from different rings may vary because of the approximate nature of the finite element 
solution. Strong variation in these estimates, commonly called domain dependence or contour 
dependence, indicates a need for mesh refinement (provided that the problem is suitable for 
contour integrals). Numerical tests suggest that the estimate from the first ring of elements 
abutting the crack front does not provide a high accuracy result so at least two contours are 
recommended. In the present study, five contours were calculated and the average was taken as 
the final ./-integral value. The method is quite robust in the sense that accurate contour integral 
estimates are usually obtained even with quite coarse meshes. 

Sharp cracks, where the crack faces lie on top of one another in the undeformed 
configuration, are usually modeled using small-strain assumptions. Focused meshes for J- 
integral calculation, as shown in Figure 8, should normally be used for small-strain fracture 
mechanics evaluations. 

For linear elastic materials, the linear elastic fracture mechanics (LEFM) approach 
predicts an r °' 5 singularity near the crack tip where r is the distance from the crack tip. In finite 
element analyses, forcing the elements at the crack tip to exhibit an r~° 5 strain singularity greatly 
improves accuracy and reduces the need for a high degree of mesh refinement at the crack tip 
[26]. This r °' 5 singularity can be produced using an eight-node quadrilateral element by moving 
the mid-side nodes to the quarter-points, as noted by Barsoum [27] and Henshell and Shaw [28]. 
The needed triangular shape of the elements at the crack tip can be achieved by further 
collapsing nodes a, b, and c, as shown in Figure 9, while maintaining the r 0,5 strain singularity 
[26]. In the present study, the crack tip was modeled with a ring of collapsed quadrilateral 
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elements as shown in Figure 9. The procedure used in the finite element analysis is described as 
follows: 

1 . Collapse one side of an 8-node element so that all three nodes, a, b, and c, have the same 
geometric location (on the crack tip). 

2. Move the mid-side nodes on the sides connected to the crack tip to the 1/4 point nearest 
the crack tip. 

This procedure will create the strain singularity of r 05 so it is sufficient for linear elastic 
fracture analysis. 
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Crack front node set. 
See section A-A below. 



Figure 8. Typical focused mesh for fracture mechanics evaluation [22]. 



* -1* * a , b, r • 


Isoparamet ric space Physical space 

Figure 9. Collapsed two-dimensional element. 
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In the current study, eight-node two-dimensional plane strain elements are used. The 
crack-tip is meshed using the technique described above to evaluate the /-integral. The 
undeformed and deformed meshes for the ASTM D3165 model are shown in Figure 10. ASTM 
D3165 and ASTM D1002 configurations are followed to model geometry in the finite element 
models. Specifically, twelve elements are used through the thickness of the adherends 
corresponding to twelve plies of T300/5208 with orientation and sequence of while 

five elements are used in the finite element mesh to model the adhesive layer in the thickness 
direction. The /-integral for five different crack lengths and six different loading conditions is 
evaluated using ABAQUS™ and compared with other models described previously. 



a) Before deformation 



b) After deformation 

Figure 10. Mesh description for the ASTM D3165 model at the singularity interface (for J- 

integral calculation). 
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ASTM D1002 Stress Model 


Figure 11 shows the configuration of the ASTM D1002 specimen which is more similar 
to most field applications. The model derivations of ASTM D1002 specimen configuration are 
very similar to those for ASTM D3165 specimens. Details can be found in the literature [6]. 
The general first-order laminated plate equations, Equations (1) - (7), are applied to both the 
upper and lower adherends. The field equations for Region 1, which is outside of the overlap 
area, are 


~ I v 

e -> 

— Li 

Region 1 


Adhesive 

|SL 

t 

Adherend 

\-2' 

Region 2 


__ __ ,_P 

L 3 — 

Region 3 


Figure 11. Configuration of ASTM D 1002 specimen. 
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and the governing equations for Region 3 are 

n l x = p 

M L y = P[0(L 3 -x 3 )-w L ] 


Q : = ~P 




V 


dx 


3 J 


The oblique angle 6 in Equations (47) - (50) is: 


(46) 

(47) 

(48) 

(49) 

(50) 

(51) 
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( 52 ) 


(h u +h L )n + ri 

l x +l 2 + l 3 

The field equations for Region 2 are the same six second-order ordinary differential equations for 
ASTM D3165 as shown in Equations ( 8 ) - (13). Therefore, the stress model of the ASTM 
D 1002 specimen includes six first-order ordinary differential equations and six second-order 
ordinary differential equations with 12 variables. A total of 18 boundary conditions are obtained 
at the two ends of each region from either continuity or applied force conditions [ 6 ]. The 
symbolic solver Maple 9 is used to solve the system of equations. The strain energy release rate 
is calculated with the similar VCCT procedures as described in previous section regarding 
ASTM D3165 specimens. 


RESULTS AND DISCUSSION 

In the present study, both ASTM D3165 and ASTM D1002 were modeled analytically to 
determine the strain energy release rate using the methods described previously. The symbolic 
solver Maple 9 was utilized as the mathematical tool. The finite element models for VCCT and 
./-integral were conducted using ABAQUS™ to verify the analytical results. 

In order to demonstrate the application of the developed model, T300/5208 
(Graphite/Epoxy) with ply thickness of 0.25 mm was used for both upper and lower adherends. 
Each adherends consist of 12 plies with orientation and sequence of [ 03 / 903 ]s- The engineering 
constants of T300/5208 are E n = 181 GPa, E 22 = 10.3 GPa, Gn = 7 .17 GPa, and vu = 0.28. For 
convenience, other mechanical properties of the adherends are assumed as 

£33 = E 22 , G 13 = Gn, vi 3 = V12, v 23 = 0.35, and G 23 = £ 22 /( 2 + 2v 23 ) = 3.815 GPa 
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The adhesive is Metbond 408 with the following properties: 

E a = 0.96 GPa G a = 0.34 GPa v a = 0.4 1 

The joint dimensions of the ASTM D3165 specimen include the central overlap length 1 0 = 30 
mm, notch size l„ =1.6 mm, and adherend lengths outside the central overlap L\ = L 5 = 78.4 mm, 
and adhesive thickness ;/ = 0.2 mm. The joint dimensions of the ASTM D1002 specimen include 
an overlap length of l 0 = 30 mm and adherend lengths outside the overlap L\ = A = 80 mm, 
which includes the gripped portion of the specimen. 

When a load P = 1,000 N/m was applied, the adhesive peel and shear stresses obtained 
from both the semi-analytical model and finite element model are shown in Figures 12 and 13 for 
ASTM D3165 specimens and in Figures 14 and 15 for ASTM D1002 specimens, respectively. 
Because of symmetry, only half of the overlap region is shown in the figures. It can be seen 
from the figures that the semi-analytical solutions and finite element results are very close, 
except in the vicinity of the edge fo = 30 mm). This is because the adhesive is assumed to 
behave as an elastic material and the adhesive stress distributions are defined as being uniform 
through the thickness in the derivation. In a real joint, the adhesive would undergo plastic 
deformations and the stress at the joint ends would be reduced. 
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Figure 12. Adhesive shear stress distribution in Region 3 for the ASTM D3165 specimen. 



Figure 13. Adhesive peel stress distribution in Region 3 for the ASTM D3 165 specimen. 
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Figure 14. Adhesive shear stress distribution in Region 2 for the ASTM D1002 specimen. 



x 2 (mm) 


Figure 15. Adhesive peel stress distribution in Region 2 for the ASTM D 1002 specimen. 
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Once the stress, strain, and displacement fields were obtained, the strain energy release rates of 
the joints were estimated and the solutions from the semi-analytical models using the VCCT 
were compared with finite element models using the VCCT and finite element models using the 
./-integral. The strain energy release rates for five different crack lengths (a = 0.15 mm, 1.5 mm, 
3 mm, 6 mm, and 9 mm) under tensile load P = 200 kN/m are shown in Figures 16 and 17 for 
ASTM D3165 and D1002 specimen geometries, respectively. It can be seen that results from 
finite element models with the VCCT and ./-integral are almost identical while the developed 
semi-analytical model deviates about 10% from the finite element models for the ASTM D3165 
specimen. For the ASTM D1002 geometry, all three models provided almost identical strain 
energy release rates except for a very short pre-crack. The strain energy release rates with crack 
length a = 3 mm under various loads are shown in Figures 18 and 19 for ASTM D3165 and 
D1002 specimen geometries, respectively. As expected, the strain energy release rates from all 
three models appear to be quadratic functions of the applied load. The semi-analytical results 
correlate very well for the models with the ASTM D1002 specimen geometry. 
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Figure 16. Strain energy release rate of an ASTM D3165 specimen as a function of initial 

crack length. 
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Figure 1 7. Strain energy release rate of an ASTM D1002 specimen as a function of initial crack 

length. 
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Figure 18. Strain energy release rate of an ASTM D3165 specimen as a function of load. 
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Figure 19. Strain energy release rate of an ASTM D1002 specimen as a function of load. 
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CONCLUSION 


A semi-analytical method was developed to calculate the strain energy release rate based 
on ASTM D3165 and D1002 specimen geometries with a prescribed interface crack. The stress 
and displacement fields for the adhesively-bonded single-lap composite joint were determined 
based on laminated anisotropic plate theory. The virtual crack closure technique (VCCT) was 
applied effectively in conjunction with the analytical stress and displacement models in 
determining the strain energy release rate. Results obtained from the developed semi-analytical 
method correlated well with results from finite element models using both the VCCT and J- 
integral. Therefore, the present study has given a description of a reasonably accurate, rapid- 
solution method for calculating the strain energy release rate of an adhesively-bonded single-lap 
composite joint. Strength predictions of representative joints are anticipated by using the critical 
strain energy release rate of the adhesive/adherend interface from the developed method. 
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